Проект IV#

Создание интерактивной карты#

1. Тема карты#

🛫 Отображение выполненных работ по аэрофотосъемке с помощью Folium 🛬#

Цель:

Создание интерактивной карты, отображающую географию и некоторую статистику проведенных работ по аэрофотосъемке и воздушному лазерному сканированию за временной срез 2016-2024 годы.

Исходя из требований, создам:

  1. 4 разных тематических слоя и 2 картографических подложки

  • Cartodbpositron

  • OSM

  • слой Субъекты РФ

  • слой Работали здесь

  • слой Места съемок

  • слой Количество проектов

  1. разных способа изображения данных:

  • картограмма (choropleth)

  • катодиаграмма (маркеры с кругами)

  • кластеризация маркеров

  • маркеры

  1. другие обязательные элементы:

  • Layer Control (управление слоями)

  1. Дополнительные интерактивные инструменты:**

  • Mouse Position

  • Mini Map

  • Marker cluster

2. Изучение данных#

2.1 Импорт библиотек#

import folium
import pandas as pd
import geopandas as gpd
import osmnx as ox
import gzip as gz
import json as js


import branca.colormap as cm

from branca.element import Html, MacroElement
from folium.plugins import MousePosition
from folium.plugins import MarkerCluster
from folium.plugins import MiniMap
from folium import IFrame

2.2 Загрузка данных о проектах#

# Загрузим данные о выполненных проектах
df = gpd.read_file('2016-2024_laser_cover.shp')

# Определяем СК
df_crs = df.crs

# Проверяем СК
print("Cистема координат слоя", df_crs)
---------------------------------------------------------------------------
DataSourceError                           Traceback (most recent call last)
Cell In[2], line 2
      1 # Загрузим данные о выполненных проектах
----> 2 df = gpd.read_file('2016-2024_laser_cover.shp')
      4 # Определяем СК
      5 df_crs = df.crs

File /opt/anaconda3/lib/python3.12/site-packages/geopandas/io/file.py:294, in _read_file(filename, bbox, mask, columns, rows, engine, **kwargs)
    291             from_bytes = True
    293 if engine == "pyogrio":
--> 294     return _read_file_pyogrio(
    295         filename, bbox=bbox, mask=mask, columns=columns, rows=rows, **kwargs
    296     )
    298 elif engine == "fiona":
    299     if pd.api.types.is_file_like(filename):

File /opt/anaconda3/lib/python3.12/site-packages/geopandas/io/file.py:547, in _read_file_pyogrio(path_or_bytes, bbox, mask, rows, **kwargs)
    538     warnings.warn(
    539         "The 'include_fields' and 'ignore_fields' keywords are deprecated, and "
    540         "will be removed in a future release. You can use the 'columns' keyword "
   (...)
    543         stacklevel=3,
    544     )
    545     kwargs["columns"] = kwargs.pop("include_fields")
--> 547 return pyogrio.read_dataframe(path_or_bytes, bbox=bbox, **kwargs)

File /opt/anaconda3/lib/python3.12/site-packages/pyogrio/geopandas.py:265, in read_dataframe(path_or_buffer, layer, encoding, columns, read_geometry, force_2d, skip_features, max_features, where, bbox, mask, fids, sql, sql_dialect, fid_as_index, use_arrow, on_invalid, arrow_to_pandas_kwargs, **kwargs)
    260 if not use_arrow:
    261     # For arrow, datetimes are read as is.
    262     # For numpy IO, datetimes are read as string values to preserve timezone info
    263     # as numpy does not directly support timezones.
    264     kwargs["datetime_as_string"] = True
--> 265 result = read_func(
    266     path_or_buffer,
    267     layer=layer,
    268     encoding=encoding,
    269     columns=columns,
    270     read_geometry=read_geometry,
    271     force_2d=gdal_force_2d,
    272     skip_features=skip_features,
    273     max_features=max_features,
    274     where=where,
    275     bbox=bbox,
    276     mask=mask,
    277     fids=fids,
    278     sql=sql,
    279     sql_dialect=sql_dialect,
    280     return_fids=fid_as_index,
    281     **kwargs,
    282 )
    284 if use_arrow:
    285     meta, table = result

File /opt/anaconda3/lib/python3.12/site-packages/pyogrio/raw.py:198, in read(path_or_buffer, layer, encoding, columns, read_geometry, force_2d, skip_features, max_features, where, bbox, mask, fids, sql, sql_dialect, return_fids, datetime_as_string, **kwargs)
     59 """Read OGR data source into numpy arrays.
     60 
     61 IMPORTANT: non-linear geometry types (e.g., MultiSurface) are converted
   (...)
    194 
    195 """
    196 dataset_kwargs = _preprocess_options_key_value(kwargs) if kwargs else {}
--> 198 return ogr_read(
    199     get_vsi_path_or_buffer(path_or_buffer),
    200     layer=layer,
    201     encoding=encoding,
    202     columns=columns,
    203     read_geometry=read_geometry,
    204     force_2d=force_2d,
    205     skip_features=skip_features,
    206     max_features=max_features or 0,
    207     where=where,
    208     bbox=bbox,
    209     mask=_mask_to_wkb(mask),
    210     fids=fids,
    211     sql=sql,
    212     sql_dialect=sql_dialect,
    213     return_fids=return_fids,
    214     dataset_kwargs=dataset_kwargs,
    215     datetime_as_string=datetime_as_string,
    216 )

File pyogrio/_io.pyx:1240, in pyogrio._io.ogr_read()

File pyogrio/_io.pyx:220, in pyogrio._io.ogr_open()

DataSourceError: 2016-2024_laser_cover.shp: No such file or directory
# Посмотрим на таблицу
df.head(3)

Расшифрую наполнение таблицы:

Столбцы

значения

OBJECTID

Числовой идентификатор объекта

project_id

Внутренне название проекта

area

Площадь в км2

year

Год съемки

SK

Итоговая система координат

model

Наличие 3д модели (есть/нет)

object_id

Внутренне название объекта съемки

Shape_Leng

Не обращаем внимания

Shape_Area

Не обращаем внимания

geometry

Координаты в градусах

Для дальнейшей работы мне понадобятся центроиды, по-этому нужно спроецировать градусы в метры. Сложность в том, что данные покрывают всютерриторию России. Остановимся на Web Mercator, для наших целей точности достаточно.

# Перепроецируем в спроецированную систему координат, выберем Web Mercator
df_webmer = df.to_crs(epsg=3857)

# Создаем GeoDataFrame с центроидами
centroids = df_webmer.copy()
centroids['geometry'] = centroids.geometry.centroid  # заменяем геометрию на центроиды

# Отображаем полигоны и центроиды на одной карте
m = df_webmer.explore(color='red', name='Полигоны съемки')  # полигоны
centroids.explore(m=m, 
                  color='blue', 
                  marker_kwds={'radius': 2}, 
                  name='Центроиды'
                  )  # точки

# Покажем слои
folium.LayerControl().add_to(m)
m
centroids.head(3)
# Переведем в систему координат обратно WGS84
centroids_wgs84 = centroids.to_crs(epsg=4326)

2.3 Загрузка данных о субъектах РФ#

Так как изначально файл превышает 25 МБ, он был заархивирован и, затем, преобразован в GeoDataFrame

# Путь к сжатому файлу
input_path = "russia_region.geojson.gz"

# Чтение и преобразование в GeoDataFrame
with gz.open(input_path, "rt", encoding="utf-8") as f: # "rt" - режим чтения текста (r = read, t = text)
    geojson_data = js.load(f)

# Преобразуем в GeoDataFrame
admin = gpd.GeoDataFrame.from_features(geojson_data["features"])

# Указываем систему координат, если она известна (например, WGS84)
admin.set_crs("EPSG:4326", inplace=True)  # ← Замените, если ваша СК другая

# Проверки
print("Система координат слоя:", admin.crs)
print("Тип данных слоя:", admin.geom_type.unique())
# Отобразим данные на карте

# Создадим базовую карту
m = folium.Map(tiles='CartoDB positron')

# Добавим слой на карту

admin.explore(m=m, color='lightblue', opacity=0.8, name='субъекты РФ')  # полигоны
df.explore(m=m, color='red', name='данные съемки')
centroids_wgs84.explore(m=m, color='blue', marker_kwds={'radius': 2}, name='Центроиды')  # точки

# Добавляем контроль слоев для возможности включения/выключения
folium.LayerControl().add_to(m)

# Автоматически подстраиваем масштаб под все объекты
m.fit_bounds(m.get_bounds())

# Отображаем карту
m

3. Визуализация данных#

Мы создадим хороплетную карту, где будет показано сколько проектов снято в регионе РФ. Для этого нам понадобится сделать пространственный join центроидов и регионов.

# Делаем пространственный join: для каждой точки подбираем тот регион , в который она попадает
#    predicate='within' означает, что точка должна лежать внутри полигона округа
join_map = gpd.sjoin(centroids_wgs84, admin, how='left', predicate='within')

# Переименовываем столбцы, чтобы с ними было удобнее работать 
join_map = join_map.rename(columns={
    'name': 'region_name',        # новое имя для столбца "название округа"
})

join_map.head(5)
join_map.columns
project_counts = join_map.groupby('region_name').agg(
    project_count=('region_name', 'size'),
    total_area=('area', 'sum')  # округляем до двух знаков после запятой
).reset_index()
project_counts['total_area'] = project_counts['total_area'].round(0).astype(int)  # округляем до целого числа
project_counts
# # Группировка по region_name и подсчет количества проектов
# project_counts = join_map.groupby('region_name').size().reset_index(name='project_count')
# project_counts

Так как при пространственном объединении (sjoin), геометрия берется от левой таблицы (centroids_wgs84), то есть — точки и join_map содержит геометрию точек, а не полигонов, нужно снова присоединить информацию о полигонах (регионов) к аггрегированной таблице project_counts чтобы построить хороплетную карту.

# Проверяем тип данных
print("Тип данных слоя", join_map.geom_type.unique())
# Присоединяем геометрию регионов обратно (из admin)
choropleth_data = admin.merge(project_counts, left_on='name', right_on='region_name', how='left')
# Проверяем тип данных
print("Тип данных слоя", choropleth_data.geom_type.unique())
# Если у каких-то регионов нет точек в project_count будет NaN.
# Заменим нулями и переведем в целочисленный тип
choropleth_data['project_count'] = choropleth_data['project_count'].fillna(0).astype(int)
choropleth_data['project_count'].isna().sum()
# Удалим лишние столбцы
choropleth_data = choropleth_data.drop(columns=['NUMPOINTS', 'region_name'])
# Заменим NaN значения в столбце total_area на 0
choropleth_data['total_area'] = choropleth_data['total_area'].fillna(0).astype(int)
choropleth_data.head()

Отметим, что для визуализации будем использовать name, потому что там полные данные

3.1 Создаем хороплетную карту#

# Создаем карту с пустым фоном
m = folium.Map(location=[68.13, 98.43], 
               zoom_start=2.5, 
               tiles=None )

# Добавляем подложки как отключаемые слои
folium.TileLayer(
    tiles='cartodbpositron',
    name='Cartodbpositron',
    control=True
).add_to(m)

folium.TileLayer(
    tiles='OpenStreetMap',
    name='OpenStreetMap',
    control=True
).add_to(m)

# Добавляем слой с субъектами
folium.GeoJson(
    choropleth_data,
    name='Субьекты РФ',
    style_function=lambda feature: {
        'fillColor': '#4CAFE5',
        'color': '#FAFAFA',
        'weight': 1.5,
        'fillOpacity': 0.4
    },
    tooltip=folium.GeoJsonTooltip(
        fields=['name'],
        aliases=[''],
        localize=True,
        sticky=True,
        style="padding-left: 4px; text-align: left;"
    )
).add_to(m)

# Добавляем слой: работали здесь (только > 0)
folium.GeoJson(
    choropleth_data[choropleth_data['project_count'] > 0],
    name='Работали здесь',
    style_function=lambda feature: {
        'fillColor': '#013A63',
        'color': '#FAFAFA',
        'weight': 1.5,
        'fillOpacity': 0.7
    }
).add_to(m)
# Создаем FeatureGroup для слоя с маркерами (изначально выключен)
marker_layer = folium.FeatureGroup(name='Места съемок', show=False)

# Кластер для маркеров
marker_cluster = MarkerCluster().add_to(marker_layer)

# Добавляем точки в кластер
for idx, row in centroids_wgs84.iterrows():
    location = [row.geometry.y, row.geometry.x]
    popup_text = f"Год съемки: {row['year']}"
    
    folium.Marker(
        location=location,
        popup=popup_text,
        icon=folium.Icon(color='blue', icon='plane', prefix='fa')
    ).add_to(marker_cluster)

# Добавляем слой с маркерами на карту
marker_layer.add_to(m)
<folium.map.FeatureGroup at 0x169f3f0b0>
# Считаем центроиды (НЕ сохраняем в geometry!)
centroids = choropleth_data.geometry.centroid
choropleth_data['lon'] = centroids.x
choropleth_data['lat'] = centroids.y

# Создаем пустой FeatureGroup для кругов (по умолчанию скрыт)
circle_layer = folium.FeatureGroup(name='Количество проектов')

# Добавляем круги
for _, row in choropleth_data.iterrows():
    if row['project_count'] > 0:
        html = f"""
        <div style="text-align: left; font-weight: normal; font-size: 12px;">
            <strong>{row['name']}</strong><br>
            Проектов: {row['project_count']}<br>
            Площадь съемки: {row['total_area']} км²
        </div>
        """
        iframe = IFrame(html=html, width=200, height=70)
        popup = folium.Popup(iframe, max_width=200)

        folium.CircleMarker(
            location=[row['lat'], row['lon']],
            radius=row['project_count'] ** 0.5 * 2,
            color='#01497C',
            weight=1,
            fill=True,
            fill_color="#CD8F0B",
            fill_opacity=0.8,
            popup=popup
        ).add_to(circle_layer)

circle_layer.add_to(m)

# Добавляем мини-карту
MiniMap().add_to(m)

# Добавляем координаты и котроль слоев
MousePosition().add_to(m)

folium.LayerControl(collapsed=False).add_to(m)
# Добавим заголовок к карте
title_html = '''
     <h3 align="center" style="font-size:17px">
     <b>Выполненные работы по аэрофотосъемке и воздушному лазерному сканированию в субъектах РФ</b></h3>
     '''
m.get_root().html.add_child(folium.Element(title_html))
<branca.element.Element at 0x1253e2b10>
m
Make this Notebook Trusted to load map: File -> Trust Notebook
m.save("index.html")